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Abstract 

Properties  of  solid  energetic  materials  depend  on  a 
large  extent  on  their  crystal  structure.  Thus,  the  structure 
determines  suitability  of  a  given  compound  for  defense 
purposes.  Since  there  are  no  simple  methods  to  predict 
crystal  structures,  such  structures  become  known  only 
after  a  given  material  has  been  synthesized  and 
crystallized.  The  structures  can  be  predicted  from 
quantum  mechanical  calculations,  but  until  recently  the 
reliability  of  such  predictions  was  very  low.  This 
situation  has  changed  with  the  development  of  symmetry- 
adapted  perturbation  theory  (SAPT)  based  on  density- 
functional  theory  (DFT)  description  of  monomers,  an 
approach  known  as  SAPT (DFT).  The  SAPT (DFT) 
potentials  for  dimers  of  energetic  molecules  were  applied 
to  predictions  of  properties  of  crystals  of  such  molecules 
in  a  combined  molecular  packing,  lattice  minimization, 
and  molecular  dynamics  simulations  study.  The 
properties  of  the  cyclotrimethylene  trinitramine  (RDX) 
crystal  predicted  from  first  principles  are  in  excellent 
agreement  with  experiment  and  the  predictions  are  even 
somewhat  better  than  achieved  by  empirical  potentials 
fitted  to  the  crystal  experimental  data.  A  similar  work  on 
the  1,1  -diamino-2, 2-dinitroethylene  (FOX-7)  crystal  is  in 
progress. 

1.  Summary 

Until  very  recently,  predictions  of  a  crystal  structure 
using  only  information  about  the  atomic  arrangement  of  a 
constituent  molecule  has  been  considered  an  impossible 
task[1_4].  While  the  structure  of  various  conformations  of 
isolated  molecules  can  be  readily  determined  using 


highly-correlated  electronic  structure  methods  and  while 
effective  methods  of  optimizing  crystal  structures  are 
available,  the  remaining  component  of  such  predictions, 
the  intermolecular  potentials  determining  the  binding 
forces  within  molecular  crystals,  were  not  known 
accurately  enough.  In  principle,  these  forces  can  be 
obtained  using  electronic  structure  methods,  but  so  far 
this  would  require  inordinate  computational  resources  to 
achieve  the  requisite  accuracy.  We  demonstrated  that  a 
recently  developed  computational  method,  symmetry- 
adapted  perturbation  theory  based  on  the  density- 
functional  description  of  monomers  symmetry- adapted 
perturbation  theory  (density-functional  theory) 
[SAPT(DFT)][5],  is  sufficiently  accurate  and  numerically 
efficient  to  facilitate  this  task.  SAPT(DFT)  has  been 
shown  to  provide  as  accurate  interaction  energies  as  high- 
level  wave-function-based  methods  at  a  greatly  reduced 
computational  cost.  In  fact,  in  a  time  comparable  to 
supermolecular  DFT  calculations.  SAPT(DFT)  was  first 
used  to  compute  the  interaction  potential  for  the 
cyclotrimethylene  trinitramine  (RDX)  dimer[6].  Next, 
molecular  packing  and  lattice  energy  minimization 
methods  applying  this  potential  were  used  to  generate  a 
large  number  of  polymorphs  of  RDX  crystal  and  then 
molecular  dynamic  (MD)  simulations  were  performed  on 
the  low-energy  polymorphs.  The  lowest-energy  structure 
from  among  all  of  the  polymorphs  corresponded  to  the 
observed  crystal[7,8].  In  another  application,  SAPT(DFT)- 
based  calculations  reproduced  the  lattice  energy  of  the 
benzene  crystal  to  within  a  few  percent[7],  an 
accomplishment  thought  unachievable  using  current  ab 
initio  methods[9].  Currently,  the  SAPT(DFT)  approach  is 
being  applied  to  the  1,1 -diamino-2, 2-dinitroethylene 
(FOX-7)  crystal.  The  present  version  of  our  computer 
codes  is  applicable  also  to  molecules  larger  than  FOX-7 
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or  RDX.  For  example,  interaction  energies  for  the  dimer 
of  perylene,  containing  64  atoms,  have  been  computed 
using  SAPT(DFT).  In  the  near  future,  the  predictive 
power  of  SAPT(DFT)  and  of  SAPT(DFT)-related 
methods  will  increase  very  significantly  due  to 
methodological  developments  pursued  in  our  group. 
These  developments  include  linearly  scaling  SAPT(DFT) 
and  new  DFT  functionals.  The  linearly-scaling 
SAPT(DFT)  should  enable  calculations  for  systems 
containing  hundreds  of  atoms.  The  SAPT(DFT)  method 
should  find  important  applications  in  development  of  new 
energetic  materials,  including  crystal  design,  screening 
molecules  for  co-crystallization,  and  identification  of  low- 
energy  polymorphs. 

2.  Introduction 

Many  compounds  of  interest  to  national  defense, 
including  energetic  materials,  are  crystals  of  organic 
molecules,  bound  by  intermolecular  (van  der  Waals) 
forces.  If  these  forces  are  known,  one  can  predict 
properties  of  such  crystals  using  crystal  packing,  energy 
minimization,  and  molecular  dynamics  methods.  Since 
most  organic  molecules  are  only  weakly  polar,  the  major 
component  of  the  forces  acting  on  such  molecules  comes 
from  dispersion  interactions,  due  to  dynamic  correlation 
of  the  electronic  motions  between  two  or  more  molecules. 
All  calculations  of  correlation  components  require  use  of 
expensive  high-level  electronic  structure  methods,  too 
costly  to  be  applied  to  compounds  of  interest  in  the  field 
of  energetic  materials.  Therefore,  until  recently, 
predictions  for  such  systems  had  to  be  based  on  empirical 
potentials,  i.e.,  relatively  simple  atom-atom  functions  with 
several  parameters  fitted  to  available  experimental  data  on 
molecular  crystals.  Since  empirical  potentials  contain 
only  information  on  systems  used  for  the 
parameterization,  this  narrows  the  range  of  the  predictions 
and,  in  particular,  makes  predictions  for  notional 
materials  unreliable.  In  fact,  empirical  potentials  are  also 
applicable  only  to  systems  in  thermodynamic  states  near 
to  those  used  for  fitting  the  parameters,  so  that  for 
example  the  high-temperature  behavior  important  for 
energetic  materials  is  typically  beyond  the  range  of 
applicability.  The  predictive  power  of  theory — mainly 
based  on  empirical  potentials — was  examined  in  a  blind 
test  conducted  by  the  Cambridge  Crystallographic  Data 
Center[10]  and  found  to  be  very  low.  The  predictions  did 
not  improve  in  subsequent  tests,  in  fact,  the  success  rate 
of  the  third  test[11]  was  even  lower  than  that  of  the 
previous  ones.  Only  very  recently,  a  fourth  test  was 
performed  and  the  success  rate  has  become  somewhat 
better[12_14].  This  is  an  unsatisfactory  situation  in  view  of 
the  importance  of  such  crystals.  For  example,  energetic 
materials  are  crystals  of  large  organic  molecules  and 


polymorphism  of  drugs  is  a  major  problem  in 
pharmaceutical  industry.  This  inability  of  theory  to 
predict  properties  of  molecular  crystals  has  been 
considered  to  be  one  of  the  major  issues  in  modern 
molecular  science[1_4]. 

We  now  know  that  the  failures  of  predictions 
discussed  above  were  due  to  insufficient  accuracy  of  the 
force  fields  used.  It  is  currently  possible  to  compute 
highly  accurate  force  fields  using  ab  initio  wave- function 
(WF)  based  methods,  but  only  for  molecules  significantly 
smaller  than  those  forming  energetic  materials.  On  the 
other  hand,  although  the  DFT  approach  can  be  applied  to 
systems  containing  hundreds  of  atoms,  conventional  DFT 
methods  poorly  describe  intermolecular  interactions 
dominated  by  the  dispersion  component.  This  problem 
has  been  overcome  by  combining  the  WF  and  DFT 
methodologies  in  symmetry- adapted  perturbation  theory 
(SAPT)[15,16]  with  the  Kohn-Sham  DFT  representation  of 
monomers[5’17_22],  resulting  in  a  method  dubbed  as 
SAPT(DFT)  or  DFT-SAPT.  The  high  efficiency  of 
SAPT(DFT)  has  been  achieved 19,23_25]  by  applying  the 
density- fitting  method[26].  [For  discussions  of 
computational  scaling  of  SAPT(DFT),  see  References  23 
and  25.]  SAPT(DFT)  calculations  take  less  time  than 
supermolecular  DFT  calculations  for  up  to  about  30  atoms 
per  monomer  since  the  dimer  DFT  calculation  is  avoided. 
SAPT(DFT)  has  been  tested  on  a  broad  range  of  systems 
and  shown  to  give  results  comparable  to  those  computed 
using  high-level  WF-based  methods[5,23,27,28].  In 
particular,  SAPT(DFT)  has  been  used  to  develop 
complete  potential  surfaces  for  the  water[29]  and 
benzene[27]  dimers.  Application  of  these  surfaces  in 
calculations  of  properties  of  clusters  and  bulk  phase  gave 
results  in  excellent  agreement  with  experiments. 
SAPT(DFT)  can  be  applied  to  fairly  large  systems  such  as 
the  dimers  of  pyrene1-30-1  or  perylene[31]. 

3.  RDX  Crystal 

Our  work  on  the  RDX  crystal  culminated  in  two 
papers  published  within  the  last  year[7,8].  Although  the 
crystal  structure  of  RDX  has  been  measured[32],  we  have 
not  used  this  information  or  any  other  experimental  data. 
Thus,  our  predictions  were  made  under  essentially  the 
same  conditions  as  in  the  blind  tests.  The  only  exception 
was  that  we  used  experimental  monomer’s  geometry 
within  the  crystal[32].  However,  the  RDX  monomer 
geometry  from  ab  initio  optimizations^ 3 ]  is  very  close  to 
the  geometry  observed  in  crystals.  Thus,  the  use  of  the 
former  geometry  would  lead  to  nearly  identical  results. 

The  SAPT(DFT)  RDX  dimer  potential[6]  defines  the 
crystal  force  field  used  to  predict  the  crystal  structure. 
The  determination  of  this  structure  consisted  of  three 
steps:  molecular  packing  (program  MOLPAK[34]), 
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potential  energy  minimization  (program  WMIN[35]),  and 
molecular  dynamics  (MD)  simulations  (program  DL- 
POLY_2[36]).  The  first  step  uses  only  a  very  simple 
generic  force  field.  For  the  second  step,  we  had  to 
develop  a  simplified  version  of  the  SAPT(DFT)  potential, 
since  the  original  form  contained  terms  not  available  in 
WMIN.  This  simplified  version  was  a  sum  of  atom-atom 
potentials  of  the  standard  Buckingham-type  exp-6  form. 
In  the  last  step,  we  used  the  original  potential[6]. 

The  molecular  packing  procedure  sampled  51  most 
often  observed  coordination  geometries  in  27  crystal- 
symmetry  groups  to  produce  6,859  structures  for  each 
geometry.  For  the  500  most  dense  structures,  the  crystal 
energy  was  optimized  (within  restrictions  of  a  given  space 
group)  using  WMIN.  The  14  lo west-energy  structures 
have  been  used  as  starting  points  for  the  isothermal- 
isostress  (NsT)  MD  simulations  at  ambient  conditions. 
Crystal  symmetry  was  not  utilized  in  this  step.  For 
comparison,  analogous  calculations  were  performed  with 
the  empirical  potential  of  Sorescu,  Rice,  and  Thompson 
(SRT)[37]  (fitted  on  the  properties  of  the  RDX  crystal). 

The  last  step  of  our  procedure,  rarely  used  in  crystal 
structure  predictions,  is  very  important  for  obtaining 
accurate  lattice  energies  and  the  correct  ordering  of  the 
polymorphs,  as  illustrated  in  Figure  1.  This  figure 
compares  the  lattice  energies  for  the  14  low-energy 
conformers  predicted  by  WMIN  minimizations  and  by 
NsT-MD  simulations.  The  NsT-MD  lattice  energies  are 
all  higher  than  the  WMIN  results  (and  the  densities,  not 
shown  in  the  figure,  are  consistently  lower).  These  trends 
can  be  attributed  to  thermal  effects  introduced  in  the  MD 
simulations,  whereas  the  WMIN  calculations  represent  a 
0  K  result;  thermal  expansion  should  generate  lower 
densities  and  higher  lattice  energies. 


WMIN/SAPT(DFT)ej^ 
I _ I  NsT-MD/SAPTlDFT) 


RDX  Polymorph 

Figure  1.  Lattice  energies  (kcal/mol)  for  the  14  low-energy 
conformers  generated  by  MOLPAK/WMIN.  Black  bars  denote 
the  results  from  MOLPAK/WMIN  (OK)  calculations  using  the 
SAPT(DFT)exp-6  potential  and  the  grey  bars  denote  the  results 
from  NsT-MD  (298K,  1  atm)  simulations  using  the  SAPT(DFT) 
potential. 


Agreement  with  experiment  for  all  cell  parameters  is 
very  good,  especially  those  of  the  NsT-MD  results:  the 
deviations  from  experiment  of  the  NsT-MD  densities  and 
cell  edge  lengths  in  the  a-,  b-,  and  c-directions  are  -0.8%, 
0.4%,  0.3%,  and  0.1%,  respectively.  Also,  the  largest 
deviation  of  the  location  of  the  molecular  mass  center  is 
less  than  0.07  A.  This  near-perfect  agreement  is 
visualized  in  Figure  1  of  Reference  7.  Other  quantities 
computed  in  Reference  8  were  the  temperature  and 
pressure  dependences  of  the  RDX  crystal,  also  achieving 
a  very  good  agreement  with  experiment,  except  for  the 
anisotropy  of  the  thermal  expansion  coefficient. 

The  excellent  agreement  of  predictions  for  the  RDX 
crystal  based  on  the  SAPT(DFT)  potential  with 
experiment  suggests  that  this  method  can  be  reliably  used 
for  other  molecular  crystals.  In  fact,  very  good  agreement 
with  experimental  lattice  energy  was  also  obtained  for  the 
benzene  crystal[7].  Since  predictions  can  be  done  entirely 
from  first  principles,  this  approach  can  be  also  used  for 
notional  materials,  in  contrast  to  the  methods  using 
empirical  potentials.  We  expect  that  our  approach  will 
find  applications  in  crystal  design,  in  development  of 
novel  energetic  materials  and  (opto)  electronic  devices, 
screening  molecules  for  co-crystallization,  and 
identification  of  low-energy  polymorphs  of 
pharmaceutical  compounds. 

4.  FOX-7  Crystal 

The  1 , 1  -diamino-2, 2-dinitroethene  [(NH2)2-C=C- 
(N02)2]  (FOX-7)  molecule  is  a  fairly  recently  developed 
energetic  material  of  similar  performance  as  RDX  and  of 
low  sensitivity [38].  This  system  is  under  active 

investigations  in  defense  laboratories.  In  addition  to  the 
general  interest  in  this  system,  the  ab  initio  work  on  it  will 
be  of  importance  since  the  empirical  HCNO  force  field 
developed  by  the  Army  Research  Laboratory  (ARL) 
group,  which  is  generally  very  successful  in  predicting 
crystal  properties  of  energetic  nitroamines,  performs 
relatively  poorly  for  the  FOX-7  crystal.  The  FOX-7 
dimer  is  a  smaller  system  than  the  RDX  dimer,  however, 
it  requires  similar  scale  calculations.  The  reason  is  that 
FOX-7  crystal  has  monomers  in  two  different 
conformations.  Therefore,  instead  of  a  single  potential, 
one  has  to  develop  two  separate  potentials  to  investigate 
the  crystal  structure.  The  development  of  an  accurate 
potential  energy  surface  for  the  FOX-7  dimer  will  allow 
modeling  of  the  properties  of  FOX-7  with  much  greater 
predictability  than  it  was  so  far  possible. 

The  work  on  the  FOX-7  dimer  potential  and  crystal 
structure  simulations  is  in  progress.  The  interaction 
energies  computed  using  SAPT(DFT)  have  been  used  to 
determine  an  analytic  potential  function.  Starting  from 
the  relative  monomer  orientations  in  the  experimental  unit 
cell,  we  computed  radial  and  angular  cross  sections  of 
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more  than  1,000  FOX-7  dimer  configurations  and 
optimized  parameters  of  a  pair  potential  based  on  the 
Coulomb  terms  describing  interaction  between  partial 
charges  on  atoms  plus  Buckingham- type  exp- 6  terms. 
The  partial  charges  were  determined  by  fitting  to  ab  initio 
molecular  multipole  moments  of  the  monomer  (through 
1=6).  The  remaining  parameters  were  fitted  to 
SAPT(DFT)  interaction  energies  directly.  We  are 
currently  assessing  the  latest  version  of  the  potential  using 
NsT  simulations  of  a  FOX-7  supercell.  At  present,  the 
lattice  vectors  are  approaching  the  experimental  values, 
however,  the  “zig-zag”  structure  observed  within  the 
FOX-7  layers  does  not  appear  to  be  properly  described  by 
the  potential.  Calculations  of  more  ab  initio  points, 
corresponding  to  orientations  within  the  FOX-7  layers, 
are  underway  and  will  be  used  to  improve  the  current 
version  of  the  potential. 

5.  Linearly-Scaling  SAPT(DFT) 

With  its  scaling  as  the  fifth  power  of  system  size, 
SAPT(DFT)  can  be  applied  to  systems  with  about  60 
atoms.  For  larger  systems,  the  scaling  can  be  reduced  to 
linear  using  proper  localization  techniques.  In  the  case  of 
intermolecular  interactions,  the  localization  can  be 
achieved  in  a  specific  way,  as  proposed  in  Reference  39. 
This  method  computes  interactions  from  standard 
SAPT(DFT)  formulas  only  in  the  fairly  small  “contact 
region”,  i.e.,  where  the  atoms  of  the  two  interacting 
monomers  are  closest  to  each  other.  All  the  remaining 
atom-atom  interactions  are  computed  from  the  well- 
established  asymptotic  formulas [15,40].  The  approach  has 
already  been  developed  for  the  electrostatic  energies[39]. 
Currently  work  is  pursued  on  the  dispersion  energies, 
which  are  very  important  in  interactions  of  energetic 
molecules. 

In  an  effort  to  localize  dispersion  energy,  we  follow  a 
similar  path  to  that  of  Reference  39-sorting  the  atomic 
sites  of  monomers  into  near  (NR)  and  far-range  (FR) 
regions,  which  divides  calculation  of  dispersion  energy 
into  two  parts,  one  assigned  for  the  exact  method  and  the 
other  for  the  asymptotic  method.  To  implement  this 
division,  we  are  using  a  basis-space  partitioning  method 
developed  by  Misquitta  and  Stone[41],  where  one  applies 
constrained  density-fitting  technique,  which  differs  from 
conventional  density-fitting  techniques  used  in 
Reference  27  by  including  the  minimization  of  the 
repulsion  between  the  contributions  of  the  fitted  basis 
functions  (centered  at  different  sites)  to  the  initial 
transition  density.  This  extra  term  in  the  minimization 
suppresses  the  nonlocal  polarizabilities  to  a  large  extent. 
Furthermore,  these  unwanted  nonlocal  polarizabilites  can 
be  diminished  by  using  a  technique  developed  by  Le 
Sueur  and  Stone[42]. 
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